function [rmin,rmax] = analyticalBounds(restr,irs,phi,opt)
% Function uses analytical results from Gafarov, Meier and Montiel-Olea
% (2018) (GMM18) to find the upper and lower bounds of the identified set 
% of the IRs. This approach assumes that only one column of the rotation
% matrix Q is restricted and that the identified set is non-empty.
% In notation of GMM18, model is 
% Y_t = c + A_1*Y_t-1 + ... + A_p*Y_t-p + B*eps_t, eps_t ~ N(0,I).
% Restrictions are on columns of B rather than on rotation 
% matrix, Q (in notation of Giacomini and Kitagawa (2018)). 
% B = A0^(-1) = Sigmatr*Q.
% Inputs:
% - restr: structure representing restrictions
% - irs: structure containing IRs and possibly long-run (cumulative) IRs
% - phi: structure containing reduced-form VAR parameters
% - opt: structure containing model information and options

eqRestr = restr.eqRestr;
signRestr = restr.signRestr;
vma = irs.vma;
lrcir = irs.lrcir;
Sigmatr = phi.Sigmatr;
B = phi.B;
ivar = opt.ivar;
jshock = opt.jshock;
H = opt.H;
p = opt.p;
const = opt.const;
signNorm = opt.signNorm;

N = size(Sigmatr,1); % No. of variables in VAR
mZ = size(eqRestr,1); % No. of equality restrictions
mS = size(signRestr,1) + 1; % No. of sign restrictions (incl. normalisation)
kmax = min(mS,N - mZ - 1); % Maximum no. of active sign restrictions

Sigmatrprime = Sigmatr';
Sigma = Sigmatr*Sigmatrprime;
Sigmainv = Sigma\eye(N); 
iConInd = [];

if mZ > 0

    Z = zeros(N,mZ);
    
    for ii = 1:mZ

        if eqRestr(ii,4) == 1 % Restriction on A0 (i.e. B^(-1))

            Z(:,ii) = Sigmainv(eqRestr(ii,2),:)';

        elseif eqRestr(ii,4) == 2 % Restriction on A0^(-1) (i.e. B)
            
            Z(eqRestr(ii,1),ii) = 1;
            % Find indices of variables with impact IR constrained to zero.
            iConInd = [iConInd; find(eqRestr(ii,1)==ivar)];

        elseif eqRestr(ii,4) == 3 % Restriction on B^(-1)*A_l

            B = reshape(B,[N*p+const,N]);
            LagCoef = B(const+(eqRestr(ii,3)-1)*N+1:const+(eqRestr(ii,3))*N,...
                eqRestr(ii,2))';
            Z(:,ii) = (LagCoef*Sigmainv)';

        elseif eqRestr(ii,4) == 4 % Restriction on LRCIR

            Z(:,ii) = lrcir(eqRestr(ii,1),:)';

        end

    end

else

    Z = [];

end 

% Construct matrix representing sign restrictions (including normalisation).

if mS > 1

    S = zeros(N,mS);

    for ii = 1:mS-1

        if signRestr(ii,5) == 1 % Sign restriction on IR
        
            S(:,ii) = vma(signRestr(ii,1),:,signRestr(ii,3)+1)'*...
                signRestr(ii,4);
            
        elseif signRestr(ii,5) == 2 % Sign restriction on A0
            
           S(:,ii) = Sigmainv(signRestr(ii,2),:)'*signRestr(ii,4);
            
        end

    end

    if signNorm == 0

        S(:,ii+1) = Sigmainv(jshock,:)'; % Normalisation on A0

    else

        S(jshock,ii+1) = 1; % Normalisation on A0^(-1) (i.e. B)

    end

else

    if signNorm == 0

        S = Sigmainv(jshock,:)'; % Normalisation on A0

    else

        S = zeros(N,1);
        S(jshock) = 1; % Normalisation on A0^(-1) (i.e. B)

    end

end

% Compute number of combinations of active restrictions.
rNum = (mZ > 0);
kk = 1;

while kk <= kmax

    rNum = nchoosek(mS,kk) + rNum;
    kk = kk + 1;

end

% Storage
fMax = zeros(H+1,rNum,length(ivar));
fMin = fMax;

cbar = 10^8; % Set penalty term

%% k = 0 - no sign restrictions.

if mZ > 0

    r = Z; % Active restrictions (equality restrictions only)

    M = eye(N) - Sigmatrprime*r*((r'*Sigma*r)\r'*Sigmatr);
    
else
    
    M = eye(N);
    
end 

for hh = 1:H+1 % For each horizon

    for ii = 1:length(ivar) % For each variable of interest

        % Reduced-form IR of ith variable at hth horizon.
        CC = vma(ivar(ii),:,hh); 

        v = sqrt(CC*Sigmatr*M*Sigmatrprime*CC'); % Value function           

        if ~isreal(v) % If v returning as complex

            v = 0;

        end

        if v ~= 0 % If value function non-zero

            % Potential solutions.
            xPlus = Sigmatr*M*Sigmatrprime*CC'./v;
            xMinus = -xPlus;

            fMaxPlus = v - 2*(1-all(S'*xPlus >= 0))*cbar;
            fMaxMinus = -v - 2*(1-all(S'*xMinus >= 0))*cbar;

            fMinPlus = v + 2*(1-all(S'*xPlus >= 0))*cbar;
            fMinMinus = -v + 2*(1-all(S'*xMinus >= 0))*cbar;

            fMax(hh,1,ii) = max(fMaxPlus,fMaxMinus);
            fMin(hh,1,ii) = min(fMinPlus,fMinMinus);          

        elseif v == 0

            % Check whether there exists a point xstar ~= 0
            % satisfying the (active) restrictions in r and the
            % inequality restrictions not in r.

            nullr = null(r');
            nullityr = size(nullr,2);
            chk = zeros(nullityr,2);

            for jj = 1:nullityr

                xPlus = nullr(:,jj);
                xMinus = -xPlus;

                chk(jj,1) = all(S'*xPlus >= 0);
                chk(jj,2) = all(S'*xMinus >= 0);

            end

            if max(chk(:)) == 1

                fMax(hh,1,ii) = 0;
                fMin(hh,1,ii) = 0;

            else

                fMax(hh,1,ii) = -2*cbar;
                fMin(hh,1,ii) = 2*cbar;

            end

        end

    end

end


%% 0 < k <= kmax

rCount = 1; % Counter for combinations of active restrictions

for kk = 1:kmax

    % Generate vector containing column indices representing all 
    % possible combinations of k of the columns of S.
    sComb = nchoosek(1:mS,kk);

    for rr = 1:size(sComb,1) % For each set of active constraints

        rCount = rCount + 1;
        
        % Active restrictions
        r = [Z S(:,sComb(rr,:))];
        % Non-active sign restrictions.
        Sna = S(:,setdiff(1:end,sComb(rr,:)));
        
        M = eye(N) - Sigmatrprime*r*((r'*Sigma*r)\r'*Sigmatr);

        for hh = 1:H+1 % For each horizon   

            for ii = 1:length(ivar) % For each variable of interest
                
                % Reduced-form IR of ith variable at hth horizon.
                CC = vma(ivar(ii),:,hh);
                
                v = sqrt(CC*Sigmatr*M*Sigmatrprime*CC'); % Value function

                if ~isreal(v) % If v returning as complex

                    v = 0;
                    
                end
                
                if v~=0 % If value function non-zero
                
                    % Potential solutions.
                    xPlus = Sigmatr*M*Sigmatrprime*CC'./v;
                    xMinus = -xPlus;

                    fMaxPlus = v - 2*(1-all(Sna'*xPlus >= 0))*cbar;
                    fMaxMinus = -v - 2*(1-all(Sna'*xMinus >= 0))*cbar;

                    fMinPlus = v + 2*(1-all(Sna'*xPlus >= 0))*cbar;
                    fMinMinus = -v + 2*(1-all(Sna'*xMinus >= 0))*cbar;

                    fMax(hh,rCount,ii) = max(fMaxPlus,fMaxMinus);
                    fMin(hh,rCount,ii) = min(fMinPlus,fMinMinus);
                    
                elseif v == 0
                    
                    % Check whether there exists a point xstar ~= 0
                    % satisfying the (active) restrictions in r and the
                    % inequality restrictions not in r.

                    nullr = null(r');
                    nullityr = size(nullr,2);
                    chk = zeros(nullityr,2);

                    for jj = 1:nullityr

                        xPlus = nullr(:,jj);
                        xMinus = -xPlus;

                        chk(jj,1) = all(Sna'*xPlus >= 0);
                        chk(jj,2) = all(Sna'*xMinus >= 0);

                    end

                    if max(chk(:)) == 1

                        fMax(hh,rCount,ii) = 0;
                        fMin(hh,rCount,ii) = 0;

                    else

                        fMax(hh,rCount,ii) = -2*cbar;
                        fMin(hh,rCount,ii) = 2*cbar;

                    end                    
               
                end

            end

        end

    end

end

% Compute bounds as max/min over fMax/fMin, respectively, for each
% horizon and variable.
rmax = permute(max(fMax,[],2),[1 3 2]);
rmin = permute(min(fMin,[],2),[1 3 2]);

% For any variable whose IR at impact is constrained to be zero (i.e.
% variables with a zero restriction on an element of A0^(-1) / B), force
% bounds of identified set to be zero. This avoids some numerical issues
% where the value function at impact is not exactly equal to zero.
rmax(1,iConInd) = 0;
rmin(1,iConInd) = 0;

end